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Abstract 

Flows of dissipative particles driven by the peristaltic motion of a tube are numerically studied. 
A transition from a slow "unjammed" flow to a fast "jammed" flow is found through the observation 
of the flow rate at a critical width of the bottleneck of a peristaltic tube. It is also found that 
the average and fluctuation of the transition time, and the peak value of the second moment of 
the flow rate exhibit powerdaw divergence near the critical point and that these variables satisfy 
scaling relationships near the critical point. The dependence of the critical width and exponents 
on the peristaltic speed and the density is also discussed. 
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I. INTRODUCTION 



Peristalsis is a progressive wave of area contraction and expansion in a tube. Transport 
due to the peristaltic motion of a tube is one of the main transport mechanisms in biological 
systems such as the esophagus, small intestine, ureter, and so forth jl| . Peristaltic transport 
is also found in the pumping of fluids, known as peristaltic pumps, which are employed in 
medical and food engineering fields 

The study of peristaltic transport has a long history and was particularly motivated by 



the desire to understand the transport of f 



uids in the ureter. Peristaltic transport in a 



Stokes fluid was studied in the late 1960s |2|, |3| and early 1970s [i-6|. A study on the 
peristaltic transport of a micropolar fluid j3] following this direction has also been reported. 
In contrast, there have been few studies on the transport of particles under a peristaltic 



10 



llj suspended 



condition such as single particles |8|, |9[ and dilute passive solid particles 
in a Newtonian fluid. Moreover, to the best of our knowledge, the peristaltic transport of 
dense particles has never been studied, even though such a situation is commonly observed 
in flows of red blood cells, in the peristaltic pumping of corrosive sand and solid foods, and 
so forth |l|. 

Here, we consider the peristaltic transport of dry and smooth granular particles, which 
is closely related to the transport of particles through a bottleneck 



discharge of grains from a silo 



3-29 
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13| , particularly the 



29]. Through previous studies it has been established that 



there are three regimes dependin g on the linear size of the bottleneck, w. If w is sufficiently 



large, particles flow continuously [16M20l|. It is also known that the mass flow rate Q satisfies 
the Beverloo law Q oc (w — wq) b for an empirical value wo, where B is 3/2 and 5/2 for two- 



and three-dimensional systems under gravity, respectively 16 



and B is 1 for disks on a 



conveyor belt 



19 



201 ] . With decreasing size of the bottleneck, the flow 



Decomes intermittent 



owing to the formation and breakdown of an arch at the outlet 



flow stops, i.e., jamming occurs 
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24J, and finally the 



29]. Note, however, that the term jamming used in this 



context is slightly di 



in recent studies 
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'erent from that in the jamming transition of granular matter discussed 



37| . Indeed, the jamming transition can be observed even for smooth, 



i.e., frictionless, granular particles, while the jamming of a granular flow at the bottleneck 
is due to the persistent arch formation by grains, which can be observed only for rough, or 
frictional, particles. What actually corresponds to the jamming transition in the discharge 



of grains appears to be the so-called dilute-to-dense transition 2lN23|. which is observed at 
a phase boundary between continuous and intermittent regimes. 

In this paper, we report the results of our simulations on the flow rate of a peristaltic 
granular flow. We consider dry, smooth, monodisperse, and spherical particles in a three- 
dimensional sinusoidally oscillating tube. After the introduction of our numerical model in 
Sec. [Ill we present the results of our simulations in Sec. IIII1 Through the observation of the 
flow rate, we find a transition from a slow "unjammed" flow to a fast "jammed" flow occurs 
if the minimum width of the peristaltic tube is smaller than a critical value. We also find 
that the average and fluctuation of the transition time, and the peak value of the second 
moment of the flow rate exhibit power-law divergence with nontrivial exponents near the 
minimum width. Moreover, we verify the existence of scaling functions for these quantities. 
It is also found that the critical width is almost independent of the density but depends on 
peristaltic velocity. In Sec. IIV[ we discuss and summarize our results. 



II. MODEL 

We adopt the three-dimensional molecular dynamics simulation for soft-core particles 
known as the discrete element method (DEM) to simulate the peristaltic flow of spheri- 
cal granular particles. We assume that the particles are dry, smooth, and monodisperse 
spheres. This means that we can ignore tangential forces that lead to the sliding and rota- 
tion of particles as well as adhesive forces. We also assume that all the material properties 
of the particles are identical. Moreover, to idealize our setup, we ignore gravity and the 
hydrodynamic force through interstitial fluids such as air. Peristaltic motion is represented 
by the propagation of a spatially oscillating wall along the axial direction of a tube, which 
is described in detail in the following. 

The system involves N monodisperse, smooth, and mobile granular spheres with diameter 
d and mass m in a peristaltic tube with length L under a periodic boundary condition along 
the direction of propagation z of peristaltic waves. The wall of the tube consists of iV w 
embedded overlapping spheres. We ignore interactions among the particles embedded in the 
wall and the deformation of the wall. Therefore, we assume that the peristaltic wave of the 
wall is uniform and perfectly controllable. See the schematic diagram shown in Fig. [T] (a). 

We label a set of particles n p = {1, . . . , N} denoting mobile particles and II W = {iV + 
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FIG. 1. (Color online) (a) Snapshot of our simulation showing some of the parameters in our model. 
See the text for details, (b) Sequential snapshots of the simulation for w = 2.0, e = 1.83 x 1CP 1 , and 
p = 2.96 x 10" 1 . Time increases from left to right, and the peristaltic wave of the wall progresses 
upward. The system remains in the unjammed flow regime when t < r. (c) Snapshot of stationary 
state for the parameters used in (b) showing the jammed flow regime when t > r. Dark (red) 
spheres and light (green) spheres are transported particles and particles embedded in the wall, 
respectively. Wall particles facing the reader are not visualized to show the transported particles 
clearly. 

1, . . . , N + N w } denoting particles embedded in the wall, where N w is the number of particles 
embedded in the wall, with the relation II = n p U Il w . For i e U p , i.e., mobile particles, the 
equation of motion is given by 



where and fj™ are respectively the elastic and viscous contact forces acting on particle 
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(/ ij + f ij ) ) 



(1) 



i by particle j. For their explicit forms, we adopt the following model: 



(2) 
(3) 



where Q(x) is a step function with Q(x) = 1 for x > and Q(x) = otherwise, k denotes 
the spring constant, 77 is the viscosity, £^ = d — jr^j, r fj - = r { - r^-, = r^/lr^l, and 



As mentioned above, we ignore the mutual deformation of the wall. Therefore, an em- 
bedded particle % G n w in the wall can be characterized by its position in the radial direction 
r (t; z) of r i = (r^t; z i ),9 i ,z i ) in cylindrical coordinates as 



for a sinusoidal peristaltic wave of the wall with amplitude 6, wavelength A, phase velocity 
c, and average tube radius a. In contrast, 9 i and z i are set to form a nearly hexagonal lattice 
with lattice constant / on the 6-z plane. See the appendix for more precise descriptions of 9 i 
and z t . We fix / = 0.3 d throughout this paper to prevent the transported particles from pen- 
etrating the tube wall, which means that N w = 9984, although we set N = 50, 100, . . . , 600 
in our simulations. 

Hereafter, we use dimensionless quantities scaled by the units of mass m, length d, and 
time ^Jm/k. For example, the length of the tube, time, and viscosity in dimensionless forms 
are given by L/d, t^k/m, and 77/ y/mk, respectively. Note that we do not introduce any 
new notations for the dimensionless quantities in the following, and use the same notations 
as those for the dimensional quantities. 

We set a = 1.5, L = 50, A = 10, and 77 = 5.48 x 10" 3 throughout this paper. The 
controllable parameters are b, c, and N, but instead we use the minimum width of the 
peristaltic tube, i.e., the width at the bottleneck, 



v ij = v i ~ v j- We solve Eq. ([[]) numerically using the Euler method with a fixed time 



interval of 5.48 x 10 2 y/m/k for each step. 






(5) 



the strain rate 




c 



(6) 
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and the volume fraction at b = 0, 



as the control parameters in this work. Note that we do not use the volume fraction p = 
N/6(a 2 + b 2 /2)L but use p instead. This is because it is difficult to fix p in simulations and 
in real experiments when b is changed. It is necessary to change both N and b at the same 
time. Moreover, b must be changed discretely because TV is a natural number. In this study, 



the restitution coefficient of the particles, e, is given by e = exp(— ttt]/ a/2 — r/ 2 ) ~ 0.988, 
i.e., the particles are almost elastic 37|. Note, however, that such "low" inelasticity or 
dissipation is necessary to reach a steady state because the energy is input to the system 
continuously by the peristaltic motion. 



III. RESULTS 

We focus on the behavior of the flow rate J(t), whose explicit definition will be given 
in Sec. II II Al under the zero-flow-rate initial condition J(0) =0. If w is sufficiently large, 
the particles can flow without becoming stuck. We call this an unjammed flow (see Fig. [TJ 
(b)). On the other hand, if w is sufficiently small, particles become stuck at bottlenecks, 
and the stuck particles rise with the wall. We call this a jammed flow (see Fig. [TJ (c)). We 
are interested in the particle behavior for intermediate values of w. In Sec. IIII B\ we present 
the results obtained by fixing p and changing w and e. Meanwhile, the results obtained by 
changing p are given in Sec. IIII CI Finally, the results for the second moment of the flow rate 
are reported in Sec. IIII PL 

It is noteworthy that the terminology of "jammed" and "unjammed" flows is counter- 
intuitive, because a jammed flow has a larger flow rate than an unjammed flow. This can 
be understood easily by considering the process relative to a frame moving with the wall. 
Indeed, the jammed state does not exhibit any motion in this frame, while the unjammed 
state has some mobile particles moving in the backward direction. Namely, some particles 
left behind by the wall motion contribute to the negative flow rate in the unjammed phase. 
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FIG. 2. (Color online) (a) Time evolutions of the normalized flow rate Jo(t, At)/ J with J = Nc/L 
for w = 2.2448, e = 1.83 x 10" 1 , and p = 2.96 x 10 -1 . The definition of J (t,At) is given by Eq. 
([9]). The observation times At for Jo(t, At) are set to T, 2T, and 5T, where T = 5.48 is the period 
of the peristaltic wave, (b) Time evolutions of the normalized flow rates J(t)/J and Jo(t,5T)/J. 
The parameters are the same as those in (a), (c) Typical time evolution of normalized flow rate 
J I J with e = 1.83 x 1CT 1 and p = 2.96 x 1CT 1 . (d) Same plot as (c) for w = 3.24 to show the 
definition of the transition time r. See text for details. 
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A. Definition of flow rate 



The flow rate is often defined 



38 



as 



ierip 

in studies employing molecular dynamics simulations and discrete element methods, where 
Pi z {t) is the ^-component of the dimensionless momentum of particle i. In real experiments, 
however, the above definition may not be appropriate because it requires the momenta of 
all particles. Instead, the following definition is usually used: 

JAt , At) = ^wi, (9) 

where AN z '(t, At) is a signed number of particles passing a section z — z' during the time 
interval [t,t + At), with +1 counted if a particle passes from z < z' to z > z' and — 1 
counted for the opposite case. Therefore, we first check that the definition given by Eq. (jSJ) 
is in agreement with that given by Eq. ()9]). 

Figure [2] (a) plots time evolutions of the flow rate J z / =0 (t, At) normalized by J = Nc/L. 
Note that if all particles rise with the wall, i.e., v { = (0, 0, c) for every particle i, J is equal 
to J. It is clearly shown that the fluctuation of the curve decreases as At increases. More- 
over, Fig. [2] (b) shows time evolutions of the normalized flow rates J(t)/J and Jo(t, 5T)/J, 
where T = A/c = 1/e. From these figures, we confirm that Jo(t, At) converges to J(t) for 
sufficiently large At. Therefore, if one is not interested in short-time behavior, Eq. dSJ) can 
be used for data analysis. Hereafter, we adopt Eq. (JHJ) for the flow rate in the subsequent 
discussion. 



B. Transition between jammed and unjammed flow 

Figure |2] (c) shows a typical example of the time evolution of the flow rate J / J for various 
w. It is clearly shown that, while an almost stationary and slow unjammed flow continues 
for large w, a transition from a slow unjammed flow to a fast jammed flow occurs for small 
w. 

We introduce the transition time r to characterize the transition from the unjammed flow 
to the jammed flow (Fig. [2] (c)). First, we linearly fit the curve of the flow rate before a 
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transition (gentle dashed line). We again fit the curve after the transition in a similar way 
(steep dashed line). Then the transition time is defined as the time at which these two fitted 
lines cross. Note that the transition time in this definition is not the lifetime of a metastable 
state but the relaxation time from the start of the unjammed flow. Such a definition is 
commonly used to characterize the relaxation from metastable states in melting dynamics 



40 



4l|. 



39| and the so-called a-relaxation in glass transitions 

Figure [3] (a) shows the w-dependence of r for various e, which suggests the existence of a 
power law between r and w c — w. Although it may be difficult to extract a universal feature 
from Fig. [3] (a), the data may satisfy the scaling form 

r~e-tf((w c -w)/e»), (10) 

where f(x) ~ x~ a for x ~ 1 and might undergo exponential decay for x ^> 1 (see Fig. [3] 
(b)). Note that the critical width w c is determined from Eq. ( fTUj) . This is because the raw 
data obtained from the simulation, as shown in Fig. [3] (a), cannot reach the critical point 
within our simulation time. From Fig. [3] (b), we numerically obtain a — 1, v — 3/2, and 
( = 7/2 for p = 2.96 x 10" 1 . The existence of the scaling law given by Eq. ffTUl) suggests 
that this jamming transition is analogous to the conventional second-order phase transition. 

We also study the w-dependence of the fluctuation of the transition time \ T = (t 2 ) — (t) 2 . 
We again find, as shown in Figs. [3] (c) and (d), the existence of the scaling form 

Xt ~e- 2 <g((w c -w)/e»), (11) 

where g(x) ~ x~^ for x ~ 1 and might undergo exponential decay for x ^> 1. The value of 
(3 is approximately equal to 3 for p = 2.96 x 10 _1 . 

We now discuss how w c depends on e. Figure 0] (a) shows the u>-dependence of J st /J, 
where J st is the stationary flow rate of the jammed flow for w < w c or that of the unjammed 
flow for w > w c . There is a jump in each stationary flow rate J st at w — w c (e). The 
transition point w c (e) strongly depends on e as shown in Fig. H] (b). This figure shows that 
w c linearly increases with decreasing e as 

w c (e) ~ -3.75e + w max , (12) 

at least for e < 0.4. 

The jump in the stationary flow rate may imply that the jamming transition is not 
continuous but discontinuous, although the scaling forms given by Eqs. ffTUj) and ffTT|) are 
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FIG. 3. (Color online) (a) Transition time r as a function of w c — w. (b) Scaling function f(x) with 
x = (w c — w)/e u and v = 3/2; see Eq. (|10p . (c) Fluctuation of transition time \t as a function of 
w c — w. (d) Scaling function g(x) with x = (w c — w)/e u ; see Eq. (fTT]h For all figures, we use the 
density p = 2.96 x 10 _1 , and circles, squares, triangles, and inverted triangles correspond to the 
data for e = 3.65 x lO^ 1 , 1.83 x 10~ 1 ,9.13 x 10~ 2 , and 4.56 x 10" 2 , respectively. 

usually a characteristic of a continuous transition. To determine whether our jamming 
transition is continuous or discontinuous, we have investigated whether there is hysteresis in 
the vicinity of the transition point. Figure H] (c) illustrates the existence of a hysteresis loop 
of the flow rate for p = 2.96 x 10" 1 and e = 1.83 x 10" 1 , which is clear evidence that this 
jamming transition is discontinuous. In the figure, squares and circles correspond to data 
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FIG. 4. (Color online) (a) Normalized stationary flow rate J s t/J as a function of w for p = 
2.96 x 10" 1 and various e, where the dashed line is given by Eq. (I13p . The symbols are the same 
as those in Fig. [3l (b) Plot of w c as a function of e for p = 2.96 x 10 _1 . (c) Hysteresis loop of flow 
rate for p = 2.96 x lO^ 1 and e = 1.83 x 10 _1 . Circles and squares correspond to Jet/ J = and 1, 
respectively. 



obtained from simulations with the initial conditions J(0) = J and J(0) = 0, respectively. 
For J(0) = J, there is a stable and steady jammed flow for w < w' c , a transition from the 
jammed flow to the unjammed flow takes place at w = w' c , and the unjammed flow is stable 
for w > w' c . It is clear that the transition point w' c satisfies w' c > w c . 

It is well known that the transition time or relaxation time from a metastable state to 
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a stable equilibrium state diverges at a critical point for a discontinuous phase transition 
42N44j| . One commonly observes that the Vogel-Fulcher-Tammann law, or its variant 



39 



t ~ exp(— C5 7 ), is satisfied for the relaxation time, where 5 is a rescaled control parameter 
or the distance from the critical or spinodal point. Meanwhile, the power-law divergence 



t ~ 5~ a is also often observed near the spinodal point in mean-field models [39|, |42j-|44| 
and in molecular dynamics simulations for sheared colloidal systems 45[. Our system is 
another example of a system in which the transition time to a stable state exhibits power-law 
divergence at a critical point as r ~ 5~ a , with 5 = l—w/w c and a = 1 for a nonequilibrium 
discontinuous phase transition. 

We also investigate the dependence of the stationary flow rate on w in the jammed and 
unjammed flow phases. The dashed curve in Fig. [4] (a) shows the analytic relation between 
the stationary flow rate J^ uld and the minimum width in a Stokes fluid under an infinite 
wavelength limit 



J. 



fluid 



3(1 -w/a) 



J 2 + 2 2 + (I -w/a) 2 ' (13) 
where <fi = b/a = 1 — w/a. In the jammed flow phase with w < w c , the stationary flow 
rate is almost independent of w and, as expected, does not follow Eq. ([TBI . On the other 
hand, the flow rate decreases as w increases in the unjammed flow phase with w > w c . It 
is surprising that the functional form of the stationary flow rate is similar to Eq. (TT3[) for 
large e, e.g., e = 0.365 (see Fig. H] (a)), although the flow rate is larger for smaller e. This 
is counter-intuitive because the Reynolds number is usually defined as Re = a 2 c/v\ with 
kinematic viscosity v, and thus, the flow rate converges to Eq. ( TT3[) as e — > in the case of 



a Stokes flow 



11. 



46(. A means of solving this puzzle might be to consider the effect 
of compressibility. Indeed, it is known that the stationary flow rate increases as c decreases 

It is also known that a granular fluid is 



if c is greater than the sound velocity 
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IS 



compressible [49|. Of course, the existence of a role of incompressibility raises an interesting 
question: why is the behavior of a suspension flow in an incompressible fluid similar to that 
of a Stokes flow? At present, this is an open question. 



C. Density dependence 

In this subsection we investigate how the flow rate depends on the density under a fixed 
strain rate of e = 1.83 x 10 _1 . Figures [5] (a) and (b) show the normalized and unnormalized 
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FIG. 5. (Color online) (a) Normalized stationary flow rate J s t/J as a function of w for various 
p. The dashed line shows Eq. (|13p . Circles, squares, triangles, inverted triangles, and diamonds 
correspond to p = 1.48 x 10 _1 ,2.96 x 10 _1 ,4.44 x 10 _1 ,5.92 x 10 -1 , and 7.40 x 10" 1 , respectively, 
(b) Unnormalized stationary flow rate, (c) Plot of w c as a function of p. (d) Plot of a as a function 
of p. All figures are obtained for e = 1.83 x 10 _1 . 



stationary flow rates J s t/J and J st as functions of w for various p, respectively. Although 
the actual flow rate J st increases as the density increases (Fig. (b)), the normalized flow 
rate decreases with increasing density (Fig. (a)). For the jammed flow phase, w < w c , 
with relatively large p, there are some particles that do not move with the wall. Note that 
the transition point w c ~ 2.2 is almost independent of the density as shown in Fig. [5] (c). 
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FIG. 6. (Color online) (a) Snapshot of stationary state for w = 2.0, e = 1.83 x 10 -1 , and p = 
7.40 x 10 -2 . Although particles do not become stuck at bottlenecks, the stationary flow rate J s % 
reaches J. (b) Snapshot of stationary state for w = 2.0, e = 1.83 x 10" 1 , and p = 7 AO x 10 _1 . 

We also demonstrate that the exponent a is a constant for 0.15 < p < 0.60 in Fig. 
(d). For a dilute gas with p < 0.10, a is much larger than 1. This might be because the 
transport of particles is not determined by the rare collisions between particles but by the 
direct momentum transfer from the peristaltic wave of the wall (Fig. [6] (a)). On the other 
hand, a is also larger than unity for p > 0.70. This might be related to the insufficient 
free volume around the antinodes of the tube (Fig. O (b)). Particles passing through a 
bottleneck soon collide with other jammed particles around the next antinode and thus they 
are deflected. These deflected particles may affect the configuration of jammed particles 
near the bottleneck. Because we focus on the appearance of the jammed flow, we can state 
that the relationship a ~ 1 is always held for a wide range of densities. 

The dashed curve in Fig. (a) shows Eq. f)13p or the u>-dependence of the flow rate in a 
Stokes fluid. We find that the flow rate in the unjammed flow phase becomes increasingly 
close to that given by Eq. ( 1T3|) with increasing density p. This also indicates that compress- 

11 




FIG. 7. (Color online) (a) Second moment of the flow rate Xj as a function of t. We set e = 
1.8 x 1CT 1 , and curves from left to right correspond to w = 2.0,2.12,2.2,2.232,2.24, and 2.244. 
(b) Scaling function h(x) of peak values of Xji Xj ax > with x = (w c — w)/e u ; see Eq. (|14|) . Circles, 
squares, and triangles correspond to e = 1.83 x 10~ 1 ,9.13 x 10~ 2 , and 4.56 x 10~ 2 , respectively. 
Both figures are obtained for p = 2.96 x 10 _1 . 

ibility has a crucial role in causing the functional form of J st / J to deviate from Eq. (fl~3l) . It 
is also worth noting that J st / J is not constant in the jammed phase if p is sufficiently high. 



D. Second moment of flow rate 

Finally, we measure the second moment of the normalized flow rate, Xj{t) — ((<^(^) 2 ) — 
(J(t)) 2 )/J 2 . From this quantity, we might detect growing length and time scales near the 
dynamical phase transition. Indeed, Krzakala and Zdeborova have recently analyzed dy- 
namic susceptibility in order to characterize the melting dynamics leading to the equilibrium 



state of spin-glasses |43l . |44(. They demonstrated that there is a growing length scale near 
the critical point even in a discontinuous phase transition. 

Figure [7] (a) plots the time evolution of Xj(t) f° r various w. Each x,j exhibits a peak at 
an intermediate time then vanishes with increasing time. The peak value Xj ax increases as 
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w approaches w c . As shown in Fig. [7] (b), the data for x!/ ax satisfy the scaling form 



xT x -^h((w c -w)/^), (14) 
where h(x) ~ x -7 for x ~ 1 and £ = 1, 7 = 3/2. 

IV. CONCLUDING REMARKS 

We demonstrated that the flow rate is a suitable order parameter for this system and 
found that a dynamic phase transition from a slow "unjammed" flow to a fast "jammed" flow 
occurs at a critical width of the peristaltic tube. We also found that the transition times, its 
fluctuations, and the peak values of the second moments of the flow rate obey power laws as 
the minimum width approaches the critical value. Moreover, we demonstrated the existence 
of scaling functions for these quantities. It was also shown that the critical width is almost 
independent of the density but depends linearly on peristaltic velocity. Nevertheless, the 
phase transition is discontinuous and exhibits a hysteresis loop under some initial conditions. 

To the best of the authors' knowledge, this is the first report demonstrating the existence 
of a phase transition in peristaltic transport. Since the pioneering work on peristaltic trans- 
port carried out by Shapiro et al. [3]] , many works have reported how the flow rate depends 
on the amplitude ratio <p = b/a = 1 — w/2a for various fluids. In contrast, our work reveals 
that the flow rate of a granular flow exhibits a finite jump at = C = 1 — w c /2a. To 
describe this behavior using a fluid model, the authors suggest that at least compressibility 
and the local fluid-solid phase transition have to be introduced into the model. We expect 
that our finding will be confirmed by experiments in the near future. 

As mentioned in Sec. [U an analogous phase transition exists in a granular flow through 



a small bottleneck, which is called the dilute-to-dense transition 
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23|. Recently Zhong et 



al. reported [22J that the characteristic time r satisfies r ~ (w c — w)~ a with a = 1.8 ± 0.2 in 
a two-dimensional granular flow through a bottleneck. Their exponent is, of course, much 
larger than ours of a = 1. There is no reason why the same exponent should be obtained in 
different systems with different dimensions, but we expect that there is a universal law that 
explains such dynamical phase transitions. 

In this paper, we ignored many aspects of realistic granular particles such as the static 
friction among contact grains, the polydispersity of grain size, particle shape, and the effects 
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of adhesive forces and interstitial fluids. Moreover, we stress that the strain-control protocol 
used in this paper is unrealistic and that a stress-control protocol should be used for a 
realistic analysis of the peristaltic transport of grains. In this sense, our paper is only the 
first step to demonstrating the possible existence of an interesting dynamic phase transition 
for the peristaltic transport of granular particles, which might have occurred as a result of 
the oversimplified model that we used. 
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Appendix: Configuration of particles embedded in a wall 

It was explained in Sec. [Til that the peristaltic tube is constructed from particles with 
identical material properties and that their radial motion r^t; Zj) is given by Eq. (T4j). In this 
appendix, we determine the azimuth and height components (0$, of the position of each 
particle i e il w = {N + 1, . . . , N + N w }. 

In Sec. HJJ we described the configuration of the wall particles as being nearly hexagonal. 
We first explain what "nearly hexagonal" means. Figure M shows a schematic diagram of 
the positions of wall particles in the 9-z plane. Here the particles form a nearly hexagonal 
lattice with lattice constant lg in the ^-direction and lattice constant V in oblique directions. 
The configuration is perfectly hexagonal if lg = I'. In our simulation setup, however, we set 
lg 7^ I', and we consider this configuration to be nearly hexagonal. This is because of the 
periodicity in the 9- and ^-directions; if the configuration was perfectly hexagonal, the tube 
length L and average tube radius a would satisfy the relation 2ix(a + d/2)/n = L/\/3n' for 
n, n' G N. To avoid imposing this constraint, we consider the above nearly hexagonal lattice 
in this paper. 
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FIG. 8. Schematic diagram of the configuration of wall particles in the 9-z plane. Here, R = a+d/2. 



We introduce a parameter I to determine the lattice constants Iq ~ / and /' ~ /. We 
denote the number of particles in a layer (particles with identical z values) as N' w and define 
the number of layers N\ = N w /N^ as 

2n(a + d/2) 



N' 

w 

Ni = 2 



I 



V31 



(A.l) 
(A.2) 



where |_^J = max{n G N\n < x} is the integer part of x. From Eqs. ( lA.lj) and ( 1A.2I) . Iq and 
/' are given as 

2vr(a + rf/2) 



N' 

w 

, 2^ L 



(A.3) 
(A.4) 
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